graficador <- function(a, b, l, f) {
x <- seq(a, b, length.out = l)
y <- x
z <- outer(x, y, f)
z[z < 0] = 0
op <- par(bg = "white")
persp(x, y, z, theta = 30, phi = 30, expand = 0.5,
col = "lightblue", ltheta = 120, shade = 0.75,
zlab = f, xlab = formalArgs(args(f))[1],
ylab = formalArgs(args(f))[2]) -> res
}
tabla <- function(dataframe)
# Genera tablas con código de colores para un dataframe
{
gt(dataframe, rownames_to_stub = TRUE) %>%
gt_color_rows(columns = dimnames.data.frame(dataframe)[[2]],
palette = "ggsci::blue_material", domain = range(dataframe)) %>%
gt_theme_nytimes()
}
FUNCION DE PRECIO función relu
price_function <- function(p_1, p_2, k = 400)
{
return(1 - (exp(p_1/k)/(exp(p_1/k) + exp(p_2/k))))
}
graficador(3, 1000, 60, "price_function")
FUNCION DE DECISON
dicotomic_function_ps <- function(c_ps, c_xb, l_ps = 0.8) {
return((1 - c_ps) * (1 - c_xb) * 0.75 +
c_ps * (1 - c_xb) * 0.9 +
c_xb * (1 - c_ps) * 0.3 +
c_ps * c_xb * l_ps)
}
dicotomic_function_xb <- function(c_ps, c_xb, l_xb = 0.9) {
return((1 - c_ps) * (1 - c_xb) * 0.7 +
c_ps * (1 - c_xb) * 0.4 +
c_xb * (1 - c_ps) * 0.8 +
c_ps * c_xb * l_xb)
}
CONTROL DE PRECIO
price_control_ps <- function(p, m = 2000, c_ps, c_xb,
l_ps = 0.8) {
m = m * dicotomic_function_ps(c_ps, c_xb, l_ps = l_ps)
return(0 + (p <= m) * (p>=0) *((-(p - m)/m)^2))
}
price_control_xb <- function(p, m = 2000, c_ps, c_xb,
l_xb = 0.9) {
m = m * dicotomic_function_xb(c_ps, c_xb, l_xb = l_xb)
return(0 + (p <= m) * (p>=0) * (((p - m)/m)^2))
}
FUNCIONES DE VENTAS
ventas_ps <- function(p_ps, p_xb, c_ps = 1, c_xb = 0,
k = 400, m = 2000, l_ps = 0.8)
{
return(5e+07 * dicotomic_function_ps(c_ps, c_xb,
l_ps = l_ps) * price_function(p_ps, p_xb, k) *
price_control_ps(p_ps, m, c_ps, c_xb, l_ps = l_ps))
}
ventas_xb <- function(p_ps, p_xb, c_ps = 0, c_xb = 0,
k = 400, m = 2000, l_xb = 0.9)
{
return(5e+07 * dicotomic_function_xb(c_ps, c_xb,
l_xb = l_xb) * price_function(p_xb, p_ps, k) *
price_control_xb(p_xb, m, c_ps, c_xb, l_xb = l_xb))
}
FUNCIONES DE BENEFICIOS
beneficios_ps <- function(p_ps, p_xb, c_ps = 0, c_xb = 0,
k = 400, m = 2000, l_ps = 0.8)
{
ventas_ps(p_ps, p_xb, c_ps = c_ps, c_xb = c_xb, k = k,
m = m, l_ps = l_ps) * (p_ps - 100)
}
beneficios_xb <- function(p_ps, p_xb, c_ps = 0, c_xb = 0,
k = 400, m = 2000, l_xb = 0.9)
{
ventas_xb(p_ps, p_xb, c_ps = c_ps, c_xb = c_xb, k = k,
m = m, l_xb = l_xb) * (p_xb - 100)
}
graficador(1, 1000, 60, "beneficios_ps")
graficador(1, 1000, 60, "beneficios_xb")
FUNCIONES DE RESOLUCIÓN
# la siguiente función calcula el equilibrio de Nash
# del problema de Bertrand
resolucion <- function(c_ps = 0,
c_xb = 0,
k = 400,
m = 2000,
x0 = 0,
y0 = 0,
N = 100,
l_ps = 0.8,
l_xb = 0.9,
eps=10^-4)
# x0 e y0 corresponde al inicio
# N corresponde al número de iteraciones
# usaremos por defecto N = 100 por funcionar bien
{
# generamos una matrix que contenga todas las iteraciones
z = matrix(0, nrow = 2, ncol = N + 1)
z[, 1] <- c(x0, y0)# almacenamos el punto inicial
i = 0
l=TRUE#criterio de parada por diverger
aux2=TRUE #criterio de parada al converger
# iniciamos el bucle
while (all(i < N, aux2,l))
{
i = i + 1
#almacenamos el óptimo para p_p
aux = optim(z[1, i], function(x)
{
-beneficios_ps(x, z[2, i], c_ps,
c_xb, k, m, l_ps = l_ps)
})$par
#comprobamos que no sea demasiado alto
if (aux <= (m * dicotomic_function_ps(c_ps, c_xb, l_ps = l_ps))) {
z[1, i + 1] <- aux
} else{
l=FALSE
}
#almacenamos el óptimo para p_x
aux = optim(z[2, i], function(y)
{
-beneficios_xb(z[1, i + 1], y,
c_ps, c_xb, k, m, l_xb = l_xb)
})$par
if (aux <= (m * dicotomic_function_xb(c_ps, c_xb, l_xb = l_xb))) {
z[2, i + 1] <- aux
} else{
l=FALSE
}
aux2=dist(z[1, i + 1]-z[1, i])>eps
}
#nos devuelve los beneficios de cada empresa
# y los precios de las consolas del equilibrio de Nash
return(list(
'beneficios' = c(
beneficios_ps(
p_ps = z[1,
i + 1],
p_xb = z[2, i + 1],
c_ps = c_ps,
c_xb = c_xb,
k,
m,
l_ps = l_ps
),
beneficios_xb(
p_ps = z[1, i + 1],
p_xb = z[2, i + 1],
c_ps = c_ps,
c_xb = c_xb,
k,
m,
l_xb = l_xb
)
),
'precios' = z[, i + 1],
'converge' = l
))
}
SOLUCIÓN
# La siguiente función calcula la p y la q
# de un problema 2x2
solucion2 <- function(e1,e2) {
# importamos las funciones de python
np <- import("numpy", convert = FALSE)
nash <- import("nashpy", convert = FALSE)
# creamos el juego
np1<-np$array(e1)
np2<-np$array(e2)
jg<- nash$Game(np1,np2)
#obtenemos el equilibrio
equilibria <- jg$support_enumeration(non_degenerate = FALSE)
for (i in iterate(equilibria)){
}
#el último i contiene la estrategia mixta, si la hay
return(list('p' = py_to_r(i)[[1]][1],
'q' = py_to_r(i)[[2]][1]))
}
# esta función resuelve el problema completo
# devolviendonos la solución del mismo
solucion <- function(k = 400,
m = 2000,
l_ps = 0.8,
l_xb = 0.9) {
# calculamos los diferentes ejercicios de Bertrand
precios1=matrix(
nrow = 2,
ncol = 2,
dimnames = list(
c("No hace crossplay", "Hace crossplay"),
c("No hace crossplay", "Hace crossplay"))
)
precios2=precios1
beneficios1=precios1
beneficios2=precios1
for (i in 0:1) {
for (j in 0:1){
r= resolucion(
k = k,
m = m,
c_ps = j ,
c_xb = i,
l_ps = l_ps,
l_xb = l_xb
)
precios1[j+1,i+1]=r$precios[1]
precios2[j+1,i+1]=r$precios[2]
beneficios1[j+1,i+1]=round(r$beneficios[1])
beneficios2[j+1,i+1]=round(r$beneficios[2])
}}
# usamos el programa anterior para calcular
# la solución del problema 2x2
solution2x2 = solucion2(precios1, precios2)
p = solution2x2$p
q = solution2x2$q
# nos devuelve la solución junto con los
# beneficios esperados
return(
list(
'p' = p,
'q' = q,
'beneficiosesperados' = list(
'ps' = cbind(p, 1 - p) %*%
beneficios1 %*%
c(q, 1 - q),
'xb' = cbind(p, 1 - p) %*%
beneficios2 %*%
c(q, 1 - q)
),
'beneficiosps' = as.data.frame(beneficios1),
'beneficiosxb' = as.data.frame(beneficios2),
#obtenemos los precios que pone cada jugador en cada opción
'preciosps' = as.data.frame(precios1),
'preciosxb' = as.data.frame(precios2)
)
)
}
sol =solucion()
sol$beneficiosesperados
## $ps
## [,1]
## [1,] 3665411948
##
## $xb
## [,1]
## [1,] 4376168454
sol$p
## [1] 0
sol$q
## [1] 0
sol$preciosps
## No hace crossplay Hace crossplay
## No hace crossplay 420.3572 249.2645
## Hace crossplay 440.2721 439.8871
sol$preciosxb
## No hace crossplay Hace crossplay
## No hace crossplay 408.0221 413.5474
## Hace crossplay 299.0318 460.4085
sol$beneficiosps
## No hace crossplay Hace crossplay
## No hace crossplay 3063835537 460007800
## Hace crossplay 3605449303 3665411948
sol$beneficiosxb
## No hace crossplay Hace crossplay
## No hace crossplay 2747973569 2749898042
## Hace crossplay 916862215 4376168454
ANÁLISIS DE SENSIBILIDAD
# Función que realiza un mallado
# de la k y la m
sensitivity <- function(lstk, lstm)
{
#usamos las variables auxiliares
psolaux = c()# almacena las p obtenidas
qsolaux = c()# almacena las q obtenidas
pssolaux = c()# almacena los beneficios de
# ps obtenidos
xbsolaux = c()# almacena los beneficios de
# xb obtenidos
ps2solaux = c()# almacena los precios de
# ps obtenidos para c_ps=0 y c_xb=1
xb2solaux = c()# almacena los precios de
# xb obtenidos para c_ps=0 y c_xb=1
for (i in lstm)# bucle para la m
{
for (j in lstk)# bucle para la k
{ # obtenemos la solución
sol = solucion(k = j, m = i)
# almacenamos los datos que interesan
psolaux = c(psolaux, round(sol$p,
4))
qsolaux = c(qsolaux, round(sol$q,
4))
pssolaux = c(pssolaux,
round(sol$beneficiosesperados$ps,
4))
xbsolaux = c(xbsolaux,
round(sol$beneficiosesperados$xb,
4))
ps2solaux = c(ps2solaux,
round(sol$preciosps[2,
2], 2))
xb2solaux = c(xb2solaux,
round(sol$preciosxb[2,
2], 2))
}
}
# transformamos los datos a dataframes
plista = as.data.frame(matrix(psolaux,
nrow = length(lstk), ncol = length(lstm),
dimnames = list(paste("k=",
lstk), paste("m=", lstm))))
qlista = as.data.frame(matrix(qsolaux,
nrow = length(lstk), ncol = length(lstm),
dimnames = list(paste("k=",
lstk), paste("m=", lstm))))
pslista = as.data.frame(matrix(pssolaux,
nrow = length(lstk), ncol = length(lstm),
dimnames = list(paste("k=",
lstk), paste("m=", lstm))))
xblista = as.data.frame(matrix(xbsolaux,
nrow = length(lstk), ncol = length(lstm),
dimnames = list(paste("k=",
lstk), paste("m=", lstm))))
xb2lista = as.data.frame(matrix(xb2solaux,
nrow = length(lstk), ncol = length(lstm),
dimnames = list(paste("k=",
lstk), paste("m=", lstm))))
ps2lista = as.data.frame(matrix(ps2solaux,
nrow = length(lstk), ncol = length(lstm),
dimnames = list(paste("k=",
lstk), paste("m=", lstm))))
# devolvemos los diferentes dataframes
return(list('p' = plista, 'q' = qlista,
'beneficiosps' = pslista, 'beneficiosxb' = xblista,
'preciosps' = ps2lista, 'preciosxb' = xb2lista))
}
# creamos el mallado
lstk = seq(200, 600, length.out = 9)
lstm = seq(1000, 3000, length.out = 9)
# evaluamos el mallado
lista = sensitivity(lstk, lstm)
# obtenemos todas las tablas
tabla(lista[["p"]])
| m= 1000 | m= 1250 | m= 1500 | m= 1750 | m= 2000 | m= 2250 | m= 2500 | m= 2750 | m= 3000 | |
|---|---|---|---|---|---|---|---|---|---|
| k= 200 | 0 | 0 | 0 | 0 | 0.8752 | 0.848 | 0.8250 | 0.8056 | 0.7898 |
| k= 250 | 0 | 0 | 0 | 0 | 0.0000 | 0.000 | 0.8726 | 0.8519 | 0.8328 |
| k= 300 | 0 | 0 | 0 | 0 | 0.0000 | 0.000 | 0.0000 | 0.0000 | 0.8721 |
| k= 350 | 0 | 0 | 0 | 0 | 0.0000 | 0.000 | 0.0000 | 0.0000 | 0.0000 |
| k= 400 | 0 | 0 | 0 | 0 | 0.0000 | 0.000 | 0.0000 | 0.0000 | 0.0000 |
| k= 450 | 0 | 0 | 0 | 0 | 0.0000 | 0.000 | 0.0000 | 0.0000 | 0.0000 |
| k= 500 | 0 | 0 | 0 | 0 | 0.0000 | 0.000 | 0.0000 | 0.0000 | 0.0000 |
| k= 550 | 0 | 0 | 0 | 0 | 0.0000 | 0.000 | 0.0000 | 0.0000 | 0.0000 |
| k= 600 | 0 | 0 | 0 | 0 | 0.0000 | 0.000 | 0.0000 | 0.0000 | 0.0000 |
tabla(lista[["q"]])
tabla(lista[['beneficiosps']])
| m= 1000 | m= 1250 | m= 1500 | m= 1750 | m= 2000 | m= 2250 | m= 2500 | m= 2750 | m= 3000 | |
|---|---|---|---|---|---|---|---|---|---|
| k= 200 | 1507584083 | 1965591008 | 2380143772 | 2752750928 | 2652873864 | 2906595062 | 3140921257 | 3354637132 | 3551799715 |
| k= 250 | 1546131330 | 2041127306 | 2501788127 | 2926159824 | 3314630443 | 3669686131 | 3444385258 | 3695436686 | 3929660465 |
| k= 300 | 1568107579 | 2087385417 | 2579991822 | 3041689354 | 3472337098 | 3871936572 | 4241712304 | 4585411821 | 4236637708 |
| k= 350 | 1581317042 | 2116865447 | 2631766986 | 3120899396 | 3583849435 | 4018537734 | 4426469866 | 4808418765 | 5165654640 |
| k= 400 | 1589610198 | 2136159879 | 2667660319 | 3178076864 | 3665411948 | 4127766690 | 4566005972 | 4979761959 | 5370182729 |
| k= 450 | 1594821540 | 2149494742 | 2692764130 | 3218917373 | 3725420728 | 4210351901 | 4672534439 | 5113045012 | 5532264982 |
| k= 500 | 1598154405 | 2158630791 | 2710735907 | 3249232976 | 3770762571 | 4273313057 | 4755872762 | 5218676287 | 5660733657 |
| k= 550 | 1600444871 | 2164962963 | 2724059154 | 3271892288 | 3805303363 | 4322349683 | 4820988533 | 5302190893 | 5765104296 |
| k= 600 | 1601867661 | 2169741342 | 2733853579 | 3289212343 | 3832239148 | 4360797813 | 4874107560 | 5369454961 | 5849185658 |
tabla(lista[['beneficiosxb']])
| m= 1000 | m= 1250 | m= 1500 | m= 1750 | m= 2000 | m= 2250 | m= 2500 | m= 2750 | m= 3000 | |
|---|---|---|---|---|---|---|---|---|---|
| k= 200 | 1839819353 | 2361427842 | 2830308934 | 3249545230 | 2194850974 | 2395607578 | 2577488075 | 2744601586 | 2898366183 |
| k= 250 | 1901269854 | 2468318160 | 2991746669 | 3472061795 | 3909612095 | 4308629825 | 2858576629 | 3059381206 | 3243179154 |
| k= 300 | 1941150950 | 2539320374 | 3101782734 | 3626632952 | 4112964205 | 4564244623 | 4979848153 | 5365708599 | 3527160129 |
| k= 350 | 1968901134 | 2588896408 | 3179063504 | 3737922487 | 4262993529 | 4754425596 | 5215309730 | 5645356619 | 6046917677 |
| k= 400 | 1989055132 | 2625235912 | 3237404138 | 3821430326 | 4376168454 | 4900972430 | 5397094732 | 5864128281 | 6304138902 |
| k= 450 | 2004316972 | 2652984610 | 3281296274 | 3885561098 | 4464085862 | 5015532028 | 5540123887 | 6039106288 | 6512146528 |
| k= 500 | 2016297824 | 2674562585 | 3315837431 | 3936115123 | 4533714984 | 5106978440 | 5655492205 | 6180745235 | 6681047097 |
| k= 550 | 2025885318 | 2691853692 | 3343519388 | 3976700336 | 4589626178 | 5181138827 | 5749669826 | 6296532312 | 6821040583 |
| k= 600 | 2033627763 | 2706077134 | 3366207940 | 4010024095 | 4635529465 | 5242055425 | 5827721737 | 6392924016 | 6938014217 |
tabla(lista[['preciosps']])
| m= 1000 | m= 1250 | m= 1500 | m= 1750 | m= 2000 | m= 2250 | m= 2500 | m= 2750 | m= 3000 | |
|---|---|---|---|---|---|---|---|---|---|
| k= 200 | 262.54 | 290.32 | 312.87 | 331.48 | 346.95 | 359.97 | 371.05 | 380.59 | 388.87 |
| k= 250 | 274.02 | 306.94 | 334.62 | 358.20 | 378.23 | 395.51 | 410.39 | 423.46 | 434.87 |
| k= 300 | 282.32 | 319.34 | 351.29 | 379.04 | 403.09 | 424.29 | 442.75 | 459.27 | 473.89 |
| k= 350 | 288.62 | 328.83 | 364.15 | 395.57 | 423.31 | 447.84 | 469.87 | 489.49 | 507.09 |
| k= 400 | 293.51 | 336.33 | 374.72 | 409.05 | 439.89 | 467.57 | 492.61 | 515.12 | 535.52 |
| k= 450 | 297.44 | 342.45 | 383.21 | 420.17 | 453.75 | 484.18 | 511.85 | 537.17 | 560.23 |
| k= 500 | 300.68 | 347.47 | 390.30 | 429.49 | 465.45 | 498.32 | 528.41 | 556.22 | 581.63 |
| k= 550 | 303.40 | 351.69 | 396.29 | 437.39 | 475.37 | 510.47 | 542.75 | 572.72 | 600.44 |
| k= 600 | 305.64 | 355.32 | 401.43 | 444.19 | 483.92 | 520.97 | 555.28 | 587.21 | 617.02 |
tabla(lista[['preciosxb']])
| m= 1000 | m= 1250 | m= 1500 | m= 1750 | m= 2000 | m= 2250 | m= 2500 | m= 2750 | m= 3000 | |
|---|---|---|---|---|---|---|---|---|---|
| k= 200 | 273.54 | 300.44 | 322.04 | 339.78 | 354.46 | 366.76 | 377.25 | 386.25 | 394.13 |
| k= 250 | 287.86 | 320.31 | 347.22 | 369.85 | 389.01 | 405.43 | 419.60 | 431.91 | 442.73 |
| k= 300 | 298.52 | 335.58 | 367.05 | 393.95 | 417.26 | 437.50 | 455.14 | 470.86 | 484.71 |
| k= 350 | 306.74 | 347.59 | 382.89 | 413.54 | 440.64 | 464.39 | 485.49 | 504.27 | 521.04 |
| k= 400 | 313.27 | 357.20 | 395.95 | 430.06 | 460.41 | 487.34 | 511.54 | 533.22 | 552.76 |
| k= 450 | 318.52 | 365.17 | 406.73 | 443.80 | 477.12 | 507.07 | 533.99 | 558.45 | 580.76 |
| k= 500 | 322.81 | 371.77 | 415.78 | 455.54 | 491.49 | 524.07 | 553.66 | 580.73 | 605.37 |
| k= 550 | 326.50 | 377.29 | 423.56 | 465.59 | 503.91 | 538.93 | 570.82 | 600.27 | 627.40 |
| k= 600 | 329.59 | 382.16 | 430.22 | 474.35 | 514.80 | 551.95 | 586.29 | 617.64 | 646.92 |
# Función que realiza un mallado
# de la l_ps y la l_xb
sensitivity2 <- function(n,a,b)
{
lstps=rbeta(n,a,b)
#usamos las variables auxiliares
psolaux = c()# almacena las p obtenidas
qsolaux = c()# almacena las q obtenidas
pssolaux = c()# almacena los beneficios de
# ps obtenidos
xbsolaux = c()# almacena los beneficios de
# xb obtenidos
ps2solaux = c()# almacena los precios de
# ps obtenidos para c_ps=0 y c_xb=1
xb2solaux = c()# almacena los precios de
# xb obtenidos para c_ps=0 y c_xb=1
for (l_ps in lstps)# bucle para la l_ps
{ # obtenemos la solución
sol = solucion(l_ps = l_ps)
# almacenamos los datos que interesan
psolaux = c(psolaux, round(sol$p,
4))
qsolaux = c(qsolaux, round(sol$q,
4))
pssolaux = c(pssolaux,
round(sol$beneficiosesperados$ps,
4))
xbsolaux = c(xbsolaux,
round(sol$beneficiosesperados$xb,
4))
ps2solaux = c(ps2solaux,
round(sol$preciosps[2,
2], 2))
xb2solaux = c(xb2solaux,
round(sol$preciosxb[2,
2], 2))
}
# devolvemos los diferentes dataframes
return(list('p' = psolaux, 'q' = qsolaux,
'beneficiosps' = pssolaux, 'beneficiosxb' = xbsolaux,
'preciosps' = ps2solaux, 'preciosxb' = xb2solaux))
}
set.seed(2)
n=100
a=8
b=2
lista2 = sensitivity2(n,a,b)
#obtenemos todas las tablas
mean(lista2[["p"]])
## [1] 0
mean(lista2[["q"]])
## [1] 0
histo<-function(lst){
m=mean(lst)
mo=mlv(lst,method = 'shorth')
q1=quantile(lst, 0.025)
q2=quantile(lst, 0.975)
hist(lst,freq=FALSE)
abline(v =m,col = "blue",lwd = 2)
abline(v =q1,col = "red",lwd = 2)
abline(v =q2,col = "red",lwd = 2)
abline(v =mo,col = "green",lwd = 2)
return(paste(
'El cuantil para 0.025 es ' , as.character(q1),'
El cuantil para 0.975 es ',as.character(q2),'
La media es ' ,as.character(m),'
La moda es ' ,as.character(mo),'
La desviación típica es ',as.character(sd(lst))
))
}
# analizamos los beneficios de Playstation
histo(lista2[['beneficiosps']])
## [1] "El cuantil para 0.025 es 1615273846.525 \n El cuantil para 0.975 es 5198874991.925 \n La media es 3802338582.66 \n La moda es 4302440780.54 \n La desviación típica es 1006283136.00399"
# analizamos los beneficios de Xbox
histo(lista2[['beneficiosxb']])
## [1] "El cuantil para 0.025 es 3883555405.2 \n El cuantil para 0.975 es 4605210455.425 \n La media es 4374980663.61 \n La moda es 4482525143.96 \n La desviación típica es 206616743.334366"
histo(lista2[['preciosps']])
## [1] "El cuantil para 0.025 es 350.83675 \n El cuantil para 0.975 es 480.531 \n La media es 439.4743 \n La moda es 458.7956 \n La desviación típica es 37.2257921408686"
histo(lista2[['preciosxb']])
## [1] "El cuantil para 0.025 es 447.37075 \n El cuantil para 0.975 es 466.60975 \n La media es 460.3759 \n La moda es 463.4226 \n La desviación típica es 5.48209631344607"